
### Replication material for Real or Imagined?: American Urban-Rural Differences in Political Values

library(dplyr)
library(tidyr)
library(ggplot2)
library(survey)
library(srvyr)
library(broom)
library(tibble)
library(haven)
library(stargazer)
library(corrplot) 
library(purrr)
library(factoextra)
library(dotwhisker)
library(xtable)
library(questionr)
library(stringr)



############# Survey data loading and cleaning

data <- read.csv(
  here::here("values_cleaned.csv")
) %>% 
  mutate(
    INCOME = case_when(
      income == 17 ~ NA_integer_,
      TRUE ~ income
    ),
    RACE = factor(
      RACE,
      levels = c("White", "African American", "Asian American", 
                 "Hispanic", "Other/2+ race")
    ),
    EDUC = factor(
      EDUC,
      levels = c("No High School", "High School Graduate",
                 "Some College", "Completed 2 year college",
                 "Completed 4 year college", "Post Grad")
    ),
    resid_Place = factor(
      resid_Place,
      levels = c("City", "Suburb", "Small Town", "Rural Area")
    ),
    Place_Party = case_when(
      resid_Place == "City" & PARTY == "Democrat" ~ "Democrat - City",
      resid_Place == "Suburb" & PARTY == "Democrat" ~ "Democrat - Suburb",
      resid_Place == "Small Town" & PARTY == "Democrat" ~ "Democrat - Small Town",
      resid_Place == "Rural Area" & PARTY == "Democrat" ~ "Democrat - Rural Area",
      resid_Place == "City" & PARTY == "Independent" ~ "Independent - City",
      resid_Place == "Suburb" & PARTY == "Independent" ~ "Independent - Suburb",
      resid_Place == "Small Town" & PARTY == "Independent" ~ "Independent - Small Town",
      resid_Place == "Rural Area" & PARTY == "Independent" ~ "Independent - Rural Area",
      resid_Place == "City" & PARTY == "Republican" ~ "Republican - City",
      resid_Place == "Suburb" & PARTY == "Republican" ~ "Republican - Suburb",
      resid_Place == "Small Town" & PARTY == "Republican" ~ "Republican - Small Town",
      resid_Place == "Rural Area" & PARTY == "Republican" ~ "Republican - Rural Area"
    ),
    Identity_Party = case_when(
      resid_identity == "City Person" & PARTY == "Democrat" ~ "Democrat - City Person",
      resid_identity == "Suburb Person" & PARTY == "Democrat" ~ "Democrat - Suburb Person",
      resid_identity == "Small Town Person" & PARTY == "Democrat" ~ "Democrat - Small Town Person",
      resid_identity == "Country Person" & PARTY == "Democrat" ~ "Democrat - Country Person",
      resid_identity == "City Person" & PARTY == "Independent" ~ "Independent - City Person",
      resid_identity == "Suburb Person" & PARTY == "Independent" ~ "Independent - Suburb Person",
      resid_identity == "Small Town Person" & PARTY == "Independent" ~ "Independent - Small Town Person",
      resid_identity == "Country Person" & PARTY == "Independent" ~ "Independent - Country Person",
      resid_identity == "City Person" & PARTY == "Republican" ~ "Republican - City Person",
      resid_identity == "Suburb Person" & PARTY == "Republican" ~ "Republican - Suburb Person",
      resid_identity == "Small Town Person" & PARTY == "Republican" ~ "Republican - Small Town Person",
      resid_identity == "Country Person" & PARTY == "Republican" ~ "Republican - Country Person"
    ))


data <- data %>% 
  mutate(
    Place_Party = factor(
      Place_Party,
      levels = c(
        "Democrat - City", "Democrat - Suburb", "Democrat - Small Town", 
        "Democrat - Rural Area", "Independent - City", "Independent - Suburb",
        "Independent - Small Town", "Independent - Rural Area", "Republican - City",
        "Republican - Suburb", "Republican - Small Town", "Republican - Rural Area"
      )),
    Identity_Party = factor(
      Identity_Party,
      levels = c(
        "Democrat - City Person", "Democrat - Suburb Person", "Democrat - Small Town Person",
        "Democrat - Country Person", "Independent - City Person", "Independent - Suburb Person",
        "Independent - Small Town Person", "Independent - Country Person", "Republican - City Person",
        "Republican - Suburb Person", "Republican - Small Town Person", "Republican - Country Person"
      )))

data$ruralid <- ((data$ru_rural1 + data$ru_rural2)-1)/10
data$urbanid <- ((data$ru_urban1 + data$ru_urban2)-1)/10


data$ruralresent <- (-(data$resent_1+data$resent_2+data$resent_3))+30
data$ruralresent



################### RUCA codes merge
## from https://agid.acl.gov/Resources/TechDocs/RUCA_3.10_zip_code_file.xlsx

RUCA <- read.csv(here::here("2006CompleteExcelRUCAfile.csv"))

RUCA <- RUCA %>% 
  mutate(
    RUCA2 = as.factor(RUCA2.0)
  )

data$zipCode <- as.numeric(data$zipCode)
data2 <- left_join(
  data,
  RUCA,
  by = c("zipCode" = "ZIPA")
)

data2$zip <- as.numeric(data$zip)
data3 <- left_join(
  data2,
  RUCA,
  by = c("zip" = "ZIPA")
)

table(data3$RUCA2.0.x)
table(data3$RUCA2.0.y)


data <- data3 %>% 
  mutate(RUCAx = case_when(RUCA2.0.x == 1 ~ "Metro",
                           RUCA2.0.y == 1 ~ "Metro",
                           RUCA2.0.x == 1.1 ~ "Metro",
                           RUCA2.0.y == 1.1 ~ "Metro",
                           RUCA2.0.x == 2 ~ "Metro",
                           RUCA2.0.y == 2 ~ "Metro",
                           RUCA2.0.x == 2.1 ~ "Metro",
                           RUCA2.0.y == 2.1 ~ "Metro",
                           RUCA2.0.x == 3 ~ "Metro",
                           RUCA2.0.y == 3 ~ "Metro",
                           RUCA2.0.x == 4 ~ "Micro",
                           RUCA2.0.y == 4 ~ "Micro",
                           RUCA2.0.x == 4.1 ~ "Micro",
                           RUCA2.0.y == 4.1 ~ "Micro",
                           RUCA2.0.x == 4.2 ~ "Micro",
                           RUCA2.0.y == 4.2 ~ "Micro",
                           RUCA2.0.x == 5 ~ "Micro",
                           RUCA2.0.y == 5 ~ "Micro",
                           RUCA2.0.x == 5.1 ~ "Micro",
                           RUCA2.0.y == 5.1 ~ "Micro",
                           RUCA2.0.x == 5.2 ~ "Micro",
                           RUCA2.0.y == 5.2 ~ "Micro",
                           RUCA2.0.x == 6 ~ "Micro",
                           RUCA2.0.y == 6 ~ "Micro",
                           RUCA2.0.x == 6.1 ~ "Micro",
                           RUCA2.0.y == 6.1 ~ "Micro",
                           RUCA2.0.x > 6.9 ~ "Rural/Small Town",
                           RUCA2.0.y > 6.9 ~ "Rural/Small Town"
  ),
  RUCAbi = case_when(RUCA2.0.x == 1 ~ "Metro",
                     RUCA2.0.y == 1 ~ "Metro",
                     RUCA2.0.x == 1.1 ~ "Metro",
                     RUCA2.0.y == 1.1 ~ "Metro",
                     RUCA2.0.x == 2 ~ "Metro",
                     RUCA2.0.y == 2 ~ "Metro",
                     RUCA2.0.x == 2.1 ~ "Metro",
                     RUCA2.0.y == 2.1 ~ "Metro",
                     RUCA2.0.x == 3 ~ "Metro",
                     RUCA2.0.y == 3 ~ "Metro",
                     RUCA2.0.x == 4 ~ "Metro",
                     RUCA2.0.y == 4 ~ "Metro",
                     RUCA2.0.x == 4.1 ~ "Metro",
                     RUCA2.0.y == 4.1 ~ "Metro",
                     RUCA2.0.x == 4.2 ~ "Metro",
                     RUCA2.0.y == 4.2 ~ "Metro",
                     RUCA2.0.x == 5 ~ "Rural/Small Town",
                     RUCA2.0.y == 5 ~ "Rural/Small Town",
                     RUCA2.0.x == 5.1 ~ "Rural/Small Town",
                     RUCA2.0.y == 5.1 ~ "Rural/Small Town",
                     RUCA2.0.x == 5.2 ~ "Rural/Small Town",
                     RUCA2.0.y == 5.2 ~ "Rural/Small Town",
                     RUCA2.0.x == 6 ~ "Rural/Small Town",
                     RUCA2.0.y == 6 ~ "Rural/Small Town",
                     RUCA2.0.x == 6.1 ~ "Rural/Small Town",
                     RUCA2.0.y == 6.1 ~ "Rural/Small Town",
                     RUCA2.0.x > 6.9 ~ "Rural/Small Town",
                     RUCA2.0.y > 6.9 ~ "Rural/Small Town"))


## Urban-rural distribution by designation:
table(data$RUCAx)
table(data$RUCAbi)


survey <- data %>% 
  srvyr::as_survey(weight = weight)




############### Sample characteristics (weighted) - in Appendix A
###############################################################

wtd.table(x = data$RUCAbi, weights = data$weight)
wtd.table(x = data$RUCAx, weights = data$weight)
table(data$FEMALE)
wtd.table(x = data$FEMALE, weights = data$weight)
1378.264/(1378.264+1420.736)
table(data$FEMALE)
1490/(1309+1490)

table(data$EDUC)
wtd.table(x = data$EDUC, weights = data$weight)

(443.2155+395.0294)/(164.5911+466.3909+618.7647+711.0084)
table(data$EDUC)
(586+252)/(586+252+334+705+92+830)

table(data$RACE)
wtd.table(x = data$RACE, weights = data$weight)
334.7841/(334.7841+300.844+824.0345+209.3148+1130.0226)

table(data$pid7)
wtd.table(x = data$pid7, weights = data$weight)

table(data$RACE)

mean(data$Age, na.rm=TRUE)
weighted.mean(data$Age,data$weight, na.rm=TRUE)

mean(data$income)
weighted.mean(data$income,data$income, na.rm=TRUE)





################ Figure 1 #####################################
###############################################################

########## RUCA LEVELS - Metro, Micro, Small town/rural

place_values_mean <- survey %>% 
  group_by(RUCAx) %>% 
  summarise_at(vars(starts_with("PV_")), survey_mean, na.rm = TRUE) %>% 
  filter(!is.na(RUCAx)) %>% 
  reshape2::melt(by = c(RUCAx), na.rm = TRUE) %>%
  mutate(valtype = ifelse(grepl("_se", variable), "se","mean")) %>%
  mutate(variable = gsub("_se", "", variable)) %>%
  pivot_wider(
    names_from = "valtype",
    values_from = "value") %>%
  mutate(
    lwr = mean - 1.96*se,
    upr = mean + 1.96*se,
    place_values = case_when(
      variable == "PV_1" ~ "Freedom",
      variable == "PV_2" ~ "Equality",
      variable == "PV_3" ~ "Economic Security",
      variable == "PV_4" ~ "Morality",
      variable == "PV_5" ~ "Individualism",
      variable == "PV_6" ~ "Social Order",
      variable == "PV_7" ~ "Patriotism"
    ),
    place_values = factor(place_values, levels = unique(place_values)))

ggplot(
  place_values_mean, 
  aes(x = place_values, y = mean, color = RUCAx))+
  geom_pointrange(
    aes(ymin = lwr, ymax = upr),
    position = position_dodge(.7))+
  scale_color_manual(
    name = "Residence",
    values = c(
      "Metro" = "goldenrod",
      "Micro" = "yellow3",
      "Rural/Small Town" = "green4"
    )
  )+
  xlab("Values")+
  ylab("Mean")+
  labs(
    title = "Political Values",
    subtitle = "By Residence",
    caption = ""
  ) +
  scale_y_continuous(limits = c(3, 5)) +
  coord_flip()+
  guides(color = guide_legend(ncol = 2))+
  theme_bw()+
  theme(
    plot.title    = element_text(hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text(hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(size = 10, colour="black"),
    axis.title    = element_text(size = 14, colour="black", face = "bold"),
    axis.text     = element_text(size = 12, colour="black", angle = 0, hjust = 0.5),
    legend.position = "bottom"
  )


######## Corresponding Appendix C Table - Means by place
export <- place_values_mean %>% 
  select(place_values, mean, se, lwr, upr)
print(xtable(export), include.rownames = FALSE)




############### FIGURE 2 - Place identity --> values ##############
###################################################################

###### Regressions ######################

v1pat <- svyglm(PV_1 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v1pat)

v1pat <- svyglm(PV_1 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v1pat)

v2order <- svyglm(PV_2 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v2order)

v3indiv <- svyglm(PV_3 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v3indiv)

v4moral <- svyglm(PV_4 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v4moral)

v5econ <- svyglm(PV_5 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v5econ)

v6equal <- svyglm(PV_6 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v6equal)

v7free <- svyglm(PV_7 ~ ruralid + urbanid + factor(RUCAx) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v7free)


plot1 <- dwplot(list(v1pat, v2order, v3indiv, v4moral, v5econ,
                     v6equal, v7free)) + ggtitle("Place Identity Predicting Political Values, Controlling for Other Factors") +
  theme_bw(base_size=18) +
  theme(panel.background = element_blank(), axis.line = element_line(colour = "black")) + xlim(-.5, 1.5) +
  scale_color_hue(labels = c('Freedom', 'Equality', 'Economic', 'Morality', 'Individualism', 'Order', 'Patriotism'), name="") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +   scale_y_discrete(breaks = c("ruralid","urbanid"), 
                                                                                   labels = c("Rural Identity", "Urban Identity"),
                                                                                   limits=c("ruralid","urbanid")) 

plot1

### Corresponding regression tables (Appendix E, Table 3)
place_values <- as.character(
  unique(place_values_mean$place_values)
) %>% str_replace_all("\n", "")

stargazer(v1pat, v2order, v3indiv, v4moral, v5econ,
          v6equal, v7free, type = "latex", header = F,
          title="Full model including place identity, residence, and controls",
          multicolumn = TRUE, column.separate = 1, font.size="small",
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))





############# Table 2 for Appendix D - Place ID with no controls
##################################################################

v1patb <- svyglm(PV_1 ~ ruralid + urbanid, design=survey)
summary(v1patb)

v2ordeb <- svyglm(PV_2 ~ ruralid + urbanid, design=survey)
summary(v2ordeb)

v3indib <- svyglm(PV_3 ~ ruralid + urbanid, design=survey)
summary(v3indib)

v4morab <- svyglm(PV_4 ~ ruralid + urbanid, design=survey)
summary(v4morab)

v5econb <- svyglm(PV_5 ~ ruralid + urbanid, design=survey)
summary(v5econb)

v6equab <- svyglm(PV_6 ~ ruralid + urbanid, design=survey)
summary(v6equab)

v7freeb <- svyglm(PV_7 ~ ruralid + urbanid, design=survey)
summary(v7freeb)

stargazer(v1patb, v2ordeb, v3indib, v4morab, v5econb,
          v6equab, v7freeb, type = "latex", header = F,
          multicolumn = TRUE, column.separate = 1,
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))



############# Table 4 for Appendix E - Place ID and controls
##################################################################


v1patx <- svyglm(PV_1 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v1patx)

v2ordex <- svyglm(PV_2 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v2ordex)

v3indix <- svyglm(PV_3 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v3indix)

v4morax <- svyglm(PV_4 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v4morax)

v5econx <- svyglm(PV_5 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v5econx)

v6equax <- svyglm(PV_6 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v6equax)

v7freex <- svyglm(PV_7 ~ ruralid + urbanid + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v7freex)

stargazer(v1patx, v2ordex, v3indix, v4morax, v5econx,
          v6equax, v7freex, type = "latex", header = F,
          multicolumn = TRUE, column.separate = 1,
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))





############# Tables 6, 7, 8 for Appendix E - split by partisanship
##################################################################

# Democrats (Table 6)
dem <- data  %>% filter(pid7 < 4) 
dems <- dem %>% 
  srvyr::as_survey(weight = weight)

v1patd <- svyglm(PV_1 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v1patd)

v2orded <- svyglm(PV_2 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v2orded)

v3indid <- svyglm(PV_3 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v3indid)

v4morad <- svyglm(PV_4 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v4morad)

v5econd <- svyglm(PV_5 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v5econd)

v6equad <- svyglm(PV_6 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v6equad)

v7freed <- svyglm(PV_7 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=dems)
summary(v7freed)

stargazer(v1patd, v2orded, v3indid, v4morad, v5econd,
          v6equad, v7freed, type = "latex", header = F,
          multicolumn = TRUE, column.separate = 1,
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))


# Republicans (Table 7)
rep <- data  %>% filter(pid7 > 4) 
reps <- rep %>% 
  srvyr::as_survey(weight = weight)

v1patr <- svyglm(PV_1 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v1patr)

v2order <- svyglm(PV_2 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v2order)

v3indir <- svyglm(PV_3 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v3indir)

v4morar <- svyglm(PV_4 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v4morar)

v5econr <- svyglm(PV_5 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v5econr)

v6equar <- svyglm(PV_6 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v6equar)

v7freer <- svyglm(PV_7 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=reps)
summary(v7freer)

stargazer(v1patr, v2order, v3indir, v4morar, v5econr,
          v6equar, v7freer, type = "latex", header = F,
          multicolumn = TRUE, column.separate = 1,
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))

# Independents (Table 8)
ind <- data  %>% filter(pid7 == 4) 
inds <- ind %>% 
  srvyr::as_survey(weight = weight)

v1pati <- svyglm(PV_1 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v1pati)

v2ordei <- svyglm(PV_2 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v2ordei)

v3indii <- svyglm(PV_3 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v3indii)

v4morai <- svyglm(PV_4 ~ ruralid + urbanid +  RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v4morai)

v5econi <- svyglm(PV_5 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v5econi)

v6equai <- svyglm(PV_6 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v6equai)

v7freei <- svyglm(PV_7 ~ ruralid + urbanid + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=inds)
summary(v7freei)

stargazer(v1pati, v2ordei, v3indii, v4morai, v5econi,
          v6equai, v7freei, type = "latex", header = F,
          multicolumn = TRUE, column.separate = 1,
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))





############# Figure 3 in main text and Table 9 for Appendix E 
############# Rural identity minus urban identity
##################################################################

data$iddiff <- data$ruralid-data$urbanid

survey <- data %>% 
  srvyr::as_survey(weight = weight)

table(data$iddiff)

## Regressions

v1patd <- svyglm(PV_1 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v1patd)

v2orded <- svyglm(PV_2 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v2orded)

v3indid <- svyglm(PV_3 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v3indid)

v4morad <- svyglm(PV_4 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v4morad)

v5econd <- svyglm(PV_5 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v5econd)

v6equad <- svyglm(PV_6 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v6equad)

v7freed <- svyglm(PV_7 ~ iddiff + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v7freed)

## Table
stargazer(v1patd, v2orded, v3indid, v4morad, v5econd,
          v6equad, v7freed, type = "latex", header = F,
          multicolumn = TRUE, column.separate = 1,
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))

## Figure 3
plotdiff <- dwplot(list(v1patd, v2orded, v3indid, v4morad, v5econd,
                        v6equad, v7freed)) + ggtitle("Identity Difference Predicting Political Values, With Controls") +
  theme(panel.background = element_blank(), axis.line = element_line(colour = "black"), text=element_text(size=16)) + xlim(-1, 2) +
  scale_color_hue(labels = c('Freedom', 'Equality', 'Economic', 'Morality', 'Individualism', 'Order', 'Patriotism'), name="") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  scale_y_discrete(breaks = c("iddiff"), 
                   labels = c("Rural ID - Urban ID"),
                   limits=c("iddiff"))
plotdiff




############# Figure and Table for Appendix F - Binary RUCA coding
##################################################################

## Figure

place_values_mean2 <- survey %>% 
  group_by(RUCAbi) %>% 
  summarise_at(vars(starts_with("PV_")), survey_mean, na.rm = TRUE) %>% 
  filter(!is.na(RUCAbi)) %>% 
  reshape2::melt(by = c(RUCAbi), na.rm = TRUE) %>%
  mutate(valtype = ifelse(grepl("_se", variable), "se","mean")) %>%
  mutate(variable = gsub("_se", "", variable)) %>%
  pivot_wider(
    names_from = "valtype",
    values_from = "value") %>%
  mutate(
    lwr = mean - 1.96*se,
    upr = mean + 1.96*se,
    place_values = case_when(
      variable == "PV_1" ~ "Freedom",
      variable == "PV_2" ~ "Equality",
      variable == "PV_3" ~ "Economic Security",
      variable == "PV_4" ~ "Morality",
      variable == "PV_5" ~ "Individualism",
      variable == "PV_6" ~ "Social Order",
      variable == "PV_7" ~ "Patriotism"
    ),
    place_values = factor(place_values, levels = unique(place_values)))

ggplot(
  place_values_mean2, 
  aes(x = place_values, y = mean, color = RUCAbi))+
  geom_pointrange(
    aes(ymin = lwr, ymax = upr),
    position = position_dodge(.7))+
  scale_color_manual(
    name = "Residence",
    values = c(
      "Metro" = "goldenrod",
      "Rural/Small Town" = "green4"
    )
  )+
  xlab("Values")+
  ylab("Mean")+
  labs(
    title = "Political Values",
    subtitle = "By Residence",
    caption = ""
  ) +
  scale_y_continuous(limits = c(3, 5)) +
  coord_flip()+
  guides(color = guide_legend(ncol = 2))+
  theme_bw()+
  theme(
    plot.title    = element_text(hjust = 0.5, size = 20, colour="black", face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 16, colour="black"),
    legend.title  = element_text(hjust = 0.5, size = 14, colour="black", face = "bold"),
    plot.caption  = element_text(size = 10, colour="black"),
    axis.title    = element_text(size = 14, colour="black", face = "bold"),
    axis.text     = element_text(size = 12, colour="black", angle = 0, hjust = 0.5),
    legend.position = "bottom"
  )


##Table

v1patb <- svyglm(PV_1 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v1patb)

v2ordeb <- svyglm(PV_2 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v2ordeb)

v3indib <- svyglm(PV_3 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v3indib)

v4morab <- svyglm(PV_4 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v4morab)

v5econb <- svyglm(PV_5 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v5econb)

v6equab <- svyglm(PV_6 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v6equab)

v7freeb <- svyglm(PV_7 ~ ruralid + urbanid + factor(RUCAbi) + pid7 + RACE + FEMALE + educ + INCOME + AGE + church + factor(region), design=survey)
summary(v7freeb)

place_values <- as.character(
  unique(place_values_mean$place_values)
) %>% str_replace_all("\n", "")

stargazer(v1patb, v2ordeb, v3indib, v4morab, v5econb,
          v6equab, v7freeb, type = "latex", header = F,
          title="Full model including place identity, residence, and controls",
          multicolumn = TRUE, column.separate = 1, font.size="small",
          no.space = TRUE, keep.stat = c("n", "rsq", "adj.rsq"))

